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Abstract 

We study numerically a lattice model of semiflexible homopolymers with near- 
est neighbor attraction and energetic preference for straight joints between bonded 
monomers. For this we use a new Monte Carlo algorithm, the 'Pruned-Enriched Rosen- 
bluth Method' (PERM). It is very efficient both for relatively open configurations at 
high temperatures, and for compact and frozen-in low-T states. This allows us to 
study in detail the phase diagram as a function of nn attraction e and stiffness x. It 
shows a ^-collapse line with a transition from open coils (small e) to molten compact 
globules (large e), and a freezing transition toward a state with orientational global 
order (large stiffness x). Qualitatively this is similar to a recently studied mean field 
theory (S. Doniach, T. Garel and H. Orland (1996), J. Chem. Phys. 105 (4), 1601), 
but there are important differences in details. In contrast to the mean field theory 
and to naive expectations, the 0-temperature increases with stiffness x. The freezing 
temperature increases even faster, and reaches the 0-line at a finite value of x. For 
even stiffer chains, the freezing transition takes place directly, without the formation 
of an intermediate globular state. Although being in conflict with mean field theory, 
the latter had been conjectured already by Doniach et al. on the basis of heuristic 
arguments and of low statistics Monte Carlo simulations. 

Finally, we discuss the relevance of the present model as a very crude model for 
protein folding. 
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1 Introduction 



The statistical mechanical study of protein folding |jT|] is still at its beginning. Minimal 
models try to represent its gross features by incorporating only those few ingredients that 
are supposed basic for its qualitative understanding. Mainly with this motivation, Doniach 
et al. Pi studied recently a model of semi-stiff lattice chains. In this model monomers are 
located at the sites of a simple cubic lattice. An attraction between non-bonded nearest 
neighbors was included to mimic the effect of average hydrophobicity. In order to induce an 
ordering phase transition between a random (molten) globule and a frozen configuration - 



which would mimic a uniquely folded protein an interaction was included which favored 
straight joints between bonds along the polymer backbone over rectangular joints. One way 
to interpret this, pointed out in [Q] , is to interpret each "monomer" as a a-helical turn (ca. 3 
amino acids), and to consider the ordering transition as a transition to a protein consisting 
only of a-helices. 

Depending on the temperature and the chain stiffness, three phases were indeed found 
in 0] by means of a mean field theory: an open coil at high T, a collapsed but 'molten' 
globule at intermediate T and low stiffness, and a 'frozen' state at low T and large stiffness. 
The coil-globule transition had all typical features of the ^-transition found at zero stiffness: 
it is second order (indeed, it is a tricritical phenomenon Q]), and Tq should depend only 
weakly on chain stiffness. The freezing transition aXT = Tp should however be first order. 
Although the mean field theory predicted that Tp < Tq for all values of the stiffness, it was 
conjectured in that this might indeed not be correct, and that Tp might indeed become 
larger than Tq for sufficiently stiff chains. If that is the case, then one should observe a 
direct first order transition from open coils to ordered states for very stiff chains. This was 
at least not contradicted by Monte Carlo simulations made by the same authors 0] , but the 
simulations were hardly convincing as the authors were not able to simulate sufficiently long 
and stiff chains. 

Actually, models similar to the above had been studied already much earlier ^ ^, |^, 
H, ^ |l^ as models for other semi-stiff polymers like, e.g., DNA. Baumgartner et al. and 
Mansfield ^, |^, |TU[ studied indeed melts consisting of many short semi-stiff chains. Thus 
they could not address the problem of 6'-collapse, but they showed very convincingly that 
there is a first-order ordering transition on the simple cubic lattice, while there is presumably 
a second order ordering transition in 2-d. Since the transition in 3-d was found for long chains 
and did not seem to become smoother with chain length, it is very plausible that it should 
coincide with the freezing transition found in [Q, in the limit of infinite chain length. The 
results of [D and 0, |, ^ 0] are indeed fully consistent. 

The opposite case of a single chain, but at parameters where freezing was out of range, 
was studied in [|, . These authors were mainly interested in the problem whether stiffness 
increases or decreases the theta temperature Tq. Naively one might expect that stiffness 
should make collapse less easy, and should therefore decrease Tq. A swelling with increased 
stiffness at fixed T Tg\mboxnon — stiff was indeed seen in in [|, ^ for short chains. But 
these authors were careful to point out that this might be a finite-size effect, and that the 
actual value of Tq (defined via the limit oo) might actually increase. For very long 

chains, stiffness might actually foster collapse: once a hairpin has been made, it might be 
much easier to follow a stiff chain than a completely flexible one. As we said already, the 
mean field theory of predicted that this effect should invalidate the naive argument, and 
Tg should be independent of stiffness. 

Seen as a rudimentary protein model, the above model lacks of course one essential in- 
gredient, namely heterogeneity. It is usually assumed that heterogeneity between individual 
amino acids is the main force which drives a collapsed polypeptide into a unique native 
configuration. Nevertheless, it might be that stiffness plays a similar role as heterogeneity, 
in which case the model of 0] might catch typical features of real protein folding. Indeed, 
in a recent treatment of random copolymers |1T[ the authors found a phase diagram (fig. 3 
of [Tl|) which is surprisingly similar to the one found in 0. 



To elucidate these intriguing questions, we decided to perform more extensive Monte 
Carlo simulations. A further motivation was to test a novel algorithm, the Pruned-Enriched- 
Rosenbluth-Method (PERM) developed by one of us [0. This is a chain growth algorithm 
superficially similar to the one used also in 0, but with some essential differences. It 
has proven extremely efficient in a number of problems, most of which involved however 
rather open configurations: free SAW's [O, ^-collapse of flexible chains |ff2l, and coagulation 



transition in dilute solutions [Tj], just to name a few. We wanted to see how it performs 
at very low energies and near first order phase transitions, before using it in more realistic 
studies of protein folding. 

The paper is organized as follows. After a brief description of the model and of the 
algorithm that we used (section 2), we present our numerical results concerning the transition 
to the globular phase (section 3) and the freezing transition (section 4). In section 5 we finally 
discuss these results, draw our conclusions, and point out further open problems. 



2 The Model and the Algorithm 

We represent a polymer as a self-avoiding random walk (SAW) on a simple cubic lattice 
0. Thus, monomers are placed on the lattice sites, and double occupancy of a site is 
strictly forbidden. Boundary conditions will be discussed below, as we used different ones 
for different purposes. The energy of the chain takes into account two contributions: a 
negative energy — e for each non-bonded occupied nearest-neighbor pair (this is the standard 
attraction used in simulations of 9 polymers ||15[); and a positive energy xe for each change 
of direction of the walk, i.e. for each non-collinear pair of successive bonds. The parameter 
X will be called stiffness. Depending on the interpretation of the model, it might represent 
the fact that trans bonds are energetically favored over gauche bonds, or that a helices are 
favored over random secondary structures. 

In the following we shall assume that e = fc^, or in other words we shall measure tem- 
peratures in units of e/ks where fcs is the Boltzmann constant. We will also use sometimes 
the Boltzmann factor q = e'^^^'^'^ = e^^'^ as a control parameter. 



The algorithm that we use, PERM, is described in detail in [|12|. Here we recall for com- 
pleteness its main aspects, adding some technical details which were important to simulate 
systems at very low temperature. 

The starting point of PERM is the Rosenbluth-Rosenbluth method (RR), developed 



already in 1955 |jT6|. In this method, a chain is built by adding a new monomer at each 
time step. Assume we are at step n, and the last placed monomer has free neighbors. 

If rrin-i > 1, the new monomer is placed with some probability Pn{k) in the k-th free 
neighbor site, and the algorithm continues. If not, we discard the chain and start a new 
one. Irrespective of the precise form of p„(/c) this would introduce a bias towards compact 
configurations with few free neighbors, if all generated chains were given the same weight. 
Thus each generated configuration carries a weight which compensates for this bias. In 
addition, this weight will take care of the Boltzmann weight. In the simplest case of uniform 
neighbor sampling, Pn{k) = l/m„_i, the total weight of a chain should be 

N 

Wn=1[ w^{k^) (1) 

n=l 



with 

wM = ^m„_ie-^'=/'=-^, (2) 

z being a fugacity we are free to choose. More generally, we can use any Pn{k) (provided it 
is non-zero for each allowed neighbor), and weights 

This algorithm will be optimal if we manage to keep Wn{k) constant and to avoid traps 
in which m„ = 0: in this ideal case - which is of course impossible to reach in practice - each 
attempted chain growth would succeed, and we would have perfect importance sampling. 

In real life each factor Wn{k) will fluctuate, giving roughly a lognormal distribution for 
Wn- Thus, for very long chains, the RR ensemble is dominated by rare configurations with 



very high weight, leading to serious statistical problems |]17 . 

To overcome this difficulty, PERM uses another classical idea of polymer simulations: 
enrichment [Q. Originally, this was devised as a method to overcome attrition, i.e. the 
exponential decrease of the number of successful attempts in 'simple sampling' (here, in 
contrast to the RR method, all neighbors of the last placed monomers are sampled with the 
same probability, whether they are free or not). It consists simply in replacing unsuccessful 
attempts by copies of successful ones (in this respect enrichment is similar to a genetic 
algorithm). In PERM, enrichment is implemented by monitoring the weight Wn of partially 
grown chains. If Wn exceeds some preselected upper threshold Wn , we make two or more 
copies of the chain; divide Wn by the number of copies made; place all except one onto a 
stack; and continue with the last copy. In this way the total weight Wn is exactly preserved, 
but it is more evenly spread on several configurations. This is done at every chain length n. 

The last entry to the stack is fetched if the current chain has reached its maximal length 
N, or if we "prune" it. Pruning (the opposite of enrichment) is done when the current weight 
Wn has dropped below some lower threshold W,^. If this happens, we draw a random binary 
number r„ with prob{r„ = 0} = prob{r„ = 1} = 1/2. If r„ = 1, we keep the chain but 
double its weight. If not, we discard ("prune") it, and continue with the last entry on the 
stack. If the stack is empty, we start a new chain. When the latter happens, we say we 
start a new "tour". Chains within one tour are correlated, but chains from different tours 
are strictly uncorrelated except through the dependence of Wn and W^ on previous tours. 

One can easily see that this algorithm is correct in the sense that the average partition 
sum estimate agrees exactly with the exact partition sum, provided we estimate it only 
between finished tours (i.e. at empty stack), and provided the total number of tours itself is 
uncorrelated with the partition sum. The latter would not be the case if we would stop the 
algorithm after some preset CPU time (we would not sample properly very large tours, and 
put too much weight on small ones). We should also remember that the algorithm cannot 
give unbiased estimates for free energies and for observables which are essentially based on 
free energies (such as end-to-end distances, which are derivatives of F with respect to some 
external field), if the fiuctuations of the partition sum estimates are large. This is the main 
problem at very low temperatures. 

In the code that we use, the algorithm is implemented recursively. A subroutine adds 
a monomer at a time. When we want to enrich a chain, we call the subroutine twice (or 



more, if we want to place more than one copy onto the stack), reducing accordingly the 
weight. When we want to prune a chain, we leave the subroutine without calling itself with 
probability 1/2, and else double the weight. Actually, instead of placing copies onto the 
stack at each enrichment event, we keep only one master copy which is updated at each step. 
In addition, we have an array of counters which tells us for every chain length how many 
copies of this length are still to be handled, and an array which shows the occupancy of each 
lattice site. 

This implementation is in contrast to the implementation of the enrichment idea in 0, 
which is more in spirit of a genetic algorithm: at each time, a large population of copies 
is kept in memory, each of the same length n, and pruning and enrichment are done by 
replacing low-weight chains by copies of high-weight chains. This can be efficient on large 
parallel machines, but it poses formidable storage and data transfer problems. For finite-size 
populations there are also corrections which are not easy to analyze. 

A crucial advantage of PERM is the fact that all controls - the selection probability 
Pn{k), the thresholds and W^, and the number of copies made - are dummies which 
can be modified at any stage of the run. We can thus adjust them automatically in response 
to problems which might arise. This is a major improvement over other recursive chain 
growth methods [|l^, |l9l , where the fugacity is a control parameter which cannot be changed 
during the run without introducing a bias. 

A first good choice is to take and proportional to the current average of Wn, 
W> = c+{Wn) and W< = c_(M4), with coefficients c± = 0(1) [0. Notice that with this 
choice the dependence on the fugacity z drops out, and we can choose it arbitrarily. For 
temperatures close to Tg we used this with c+ = 3 to 10, and c^/c_ ~ 10 to 50. Within this 
wide range of parameters this was very efficient. 

Problems arose however at very low T, since there the current estimate (Wn) might be 
very wrong. The most dangerous and common situation is that we underestimate (Wn) since 
we have not yet encountered a "good" (i.e., low energy) configuration. When such a configu- 
ration finally appears, our thresholds are too low (occasionally by many orders of magnitude), 
and we produce an enormous amount of copies without pruning them sufficiently. Our way 
out of this dilemma consists in letting c± depend on the total number of copies of length 
n already made during the present tour. If this number is becoming too large, we reduce 
c±. But this alone would reduce the total number of very long chains produced, since these 
are most effected by fluctuations. The algorithm is most efficient if the sample size is the 
same for all chain lengths. We cannot enforce this precisely, but we can place a strong bias 
towards it by choosing c± proportional to the total sample size. For technical reasons, we 
indeed segmented the set of tours into bunches of typical size 10^ to 10^. We replaced the 
number of conflgurations made during the present tour by that made during the present 
bunch. Let us call this M^""'^'^, while M*"**^' is the total number of conflgurations generated 
so far. We then used 

J^) (l + Mr-Va)' (4) 

with c^^ 5, c*^/c° =~ 20, a 10^ — 10^, and exponents a and (3 either 1 or 2. The latter 
was needed for the lowest values of T. 

The prefactors c\ can be learned from preliminary runs with small chain lengths. A more 
systematic strategy which was quite successful consisted in a rudimentary genetic algorithm 



(with mutation and replacement by the fittest only, but without cross-over) by which we let 
a population of pairs (c° , c° ) evolve. 

At high temperatures, it was sufficient to make just one new copy per enrichment event. 
At low T this was not enough to prevent weights from growing too large. We used a number 
of schemes which all allowed a large number of copies, but only if Wn grew excessively large. 



A good choice was ricopy = int[Ayiy„/W>]. 

The last point which can change the efficiency of the algorithm is the neighbor selection 



probability Pn{k). In agreement with [12| we found that it was not useful to include in it 
the Boltzmann weight for the contact energy. Most likely, a large number of contacts gives 
an immediate advantage but leads to a higher risk of getting blocked later. But for large x 
we found it important to favor straight steps over right angles. In the simplest case we thus 
used 

{l/m„_i going straight is blocked 

g/32.y (^run-i + e^^ — 1) going straight 
l/(m„_i + e^^' — 1) right angle turn, although straight is not blocked 

(5) 

In the ordered phase, it was useful to go even further and enhance Pn{k) even more for 
straight steps if the previous step had also been straight. 



3 The ^-collapse line 

For small stiffness, the end-to-end distance is the easiest and most straightforward in- 
dicator for the coil-globule transition. At Tg, we expect Re^N ~ N^^'^ up to logarithmic 
corrections |TB[. But, as we noted already, this collapse is delayed for stiff chains. Thus 
Re is not very useful for large values of x. This is illustrated in fig.l, where we see that we 
need chains of length > 1000 to pin down Tq. 

A more useful observable could be the second virial coefficient which should vanish at 



Tg. But measuring it for interacting chains is a bit awkward |]15[. For a-thermal chains there 
exists the very elegant Karp-Luby algorithm [^, but no similarly elegant algorithm seems 
to be known for thermal chains. 



We found the most efficient method to measure Tg to be one due to Dickman ||2T[]. Here 
we use the fact that the pressure vanishes at Tg, and that pressure is defined as the derivative 
of the free energy with respect to volume. We thus study simultaneously two systems in 
volumes of different sizes, and compare their partition functions. More precisely, in order to 
reduce statistical fluctuations, we start with a single chain on a large lattice with periodic 
b.c. with period Li (actually we used 'helical' b.c. where each site is labeled by a single 
index i, and i + Lf = i). At each monomer insertion we check that it would still be self 
avoiding and would not change its contact energy if we would replace the periodicity by 
Lq = Li/2. As soon as this is violated, we replace the chain by two copies, one on volume 
Lf and the other on Lq, and measure the quantity 



(6) 



where is the partition sum. This vanishes exactly for short chains. For very long chains, 
it must go to —1. For T > Tg it is negative for all N (pressure is always positive), but for 
T <Tq there should exist a range where pressure is negative and thus AZ{N, Lq) > 0. For 
X — 10 this is shown in fig. 2. In this figure, the same set of temperatures are used as in 
fig.l. We see very clearly that all curves except one correspond to already collapsed chains. 
Actually, the argument is a bit more subtle due to finite size corrections which we have so 
far neglected. Simulations show that AZ{N, Lq) > even at T = Tg, but its maximum with 
respect to N does not increase with Lq unless T <Tq. Thus a precise determination of Tq is 
possible by comparing different lattice sizes. Since maxjv AZ{N, Lq) is slowly varying with 
stiffness x, it is sufficient to make this finite size analysis at few values of x. 

The results are shown in fig.3. They fully agree with those obtained from Re,N-, but are 
much more precise for large x. Most of them were obtained with Lq = 2^ and N fa 5000. 
The error on Tg is typically less than 1%, independent of x. We see clearly that Tg increases 
with stiffness, i.e. stiffness increases the tendency to collapse. 



4 The freezing transition 

At low temperature the chains are expected to undergo a first order transition toward a 
phase characterized by global order, at an x-dependent temperature Tp(x). The low energy 
configurations should appear as a bundle of linear parallel pieces in contact one to each 
other, with as few turns as possible. 

To monitor the freezing transition, and to verify that it is first order, we observed four 
different quantities: 



1. The average number of contacts per monomer, (m), which measures the contact energy. 
For large N this quantity should show a a discontinuity at Tp. But this discontinuity 
is very broad at finite chain lengths, except at very large values of x. Thus we have 
very large finite size effects, and (m) is not practically useful as an order parameter. 

2. The fraction fg of straight (trans) bonds, fs- This measures the stiffness energy and 
the local ordering of the chain. In contrast to (m), this seems very useful as order 
parameter. Above Tp it is found to be close to the naive expectation obtained for 
chains without self-avoidance and nn-attraction. 



{l + {J\f-2)q-^) ^ ^ 

(here J\f is the lattice coordination number; in our case, J\f — 6). This is a particularly 
good approximation for large x. Deviations for small x can largely be understood as 
effects of self-avoidance. In the frozen phase fs is close to 1. 

3. The fraction of bonds directed in the privileged direction. Let us denote by nx,ny, 
and Uz the bonds parallel to either of the three coordinate axes. Let us define TT-max = 
maxj=2:,y,z "^j, "^min = ^^'<^i=x,y,z ni, and p = 1 — Timin/''^max- If there is no directional 
ordering at all, we have (p) ~ 1/\/N due to the central limit theorem. In the opposite 
case of an ordered phase, we have p const for — > oo. In the intermediate case 



of weak directional ordering, a mean field type argument predicts a power law decay, 
p ~ A^~", with non-universal exponent a. 



Notice that neither Rg nor the gyration radius are very useful for detecting the ordered 
phase. We expect to see some changes at Tp since configurations should change from spherical 
globules to more rugged shapes, but we cannot expect this to be very systematic and easy 
to observe. For this reason, chain sizes were not measured for T ^ Tp, except for very stiff 
chains where Tp is rather large and where finally the collapse is without intermediate globule 
state. 

We found a dramatic dependence on x in the ease of locating Tp- In contrast to 
who were not able to go beyond x = 3, we found that the freezing transition is much easier 
to observe for large x than for small x. This should have been expected intuitively: when 
cooling down at small x, we have first a collapse to a disordered globule, and ordering sets 
in only at very low T when the mobility of the chain is extremely low. Thus we are bound 
to find important metastable states and long trapping times. For large x, in contrast, the 
ordering sets in at rather high T when the chain is still highly mobile. To illustrate this, we 
compare in figs. 4 and 5 the behavior of fs for x = 3 with that for x = 10. In both cases 
we see quite clear phase transitions, but it is much sharper for a; = 10 than for a; = 3. The 
same is true for p, see figs. 6 and 7. 

Indeed, we were not able to obtain reliable results for x = 1 and > 150, where the 
authors of claimed to see a clean ordering transition. The problem is that partition sum 
estimates fluctuate wildy in this region. At the freezing point for for a; = 1 and = 200, 
this could involve many orders of magnitude even in samples of several million chains. Thus 
even very large statistical samples were dominated by only few large-weight configurations. 
The authors of grew populations of only 20,000 chains. Our guess is that they grew 
many such populations for each set of control parameters, and averaged the results without 
weighting them with Z^- Correct averages should include this factor. Neglecting this would 
reduce greatly the statistical error estimates, at the risk of making uncontrolled systematic 
errors. 

To locate precisely a first order transition, we should in principle make a finite size scaling 
analysis. We indeed see in figs. 4 to 7 important finite size effects: with increasing A^ the 
transition point seems to shift towards higher values of T. This is as expected: for small A^ 
we have important surface effects which diminish the cooperativity of the interaction. But 
fluctuations are too large to allow a systematic analysis. 

In spite of all these problems, we were able to determine Tp in a wide range of x values 
(fig.3). Our estimates for Tp for small x agree nicely with those of and of ^ |^, |TD|. The 
ordering transition for concentrated chains (hamiltonian walks) ^, |^, |l^ should coincide 
with the presently studied transition in the limit T — > 0. From these references we expect 



Tp/x — > 0.66 for Tp —> resp. Tp/x 0.82 [|10] for Tp 0. Our data are in better 
agreement with the latter, although the large statistical errors for and the evident curvature 
of the transition line makes an extrapolation difficult. A possible alternative fit to our data 
is Tp (X x^'^^. This fit is indeed better numerically and holds for the entire range of x, but 
we see no theoretical basis for it at small values of x, and it would contradict all previous 
numerical H, ^ @] and theoretical results. 

The most conspicuous result is that Tp reaches the coil-globule transition temperature 



at X ~ 13. Beyond this triple point, we have a direct first order collapse from open coils 
to folded structures. The existence of this direct transition is also seen in R^. which drops 
suddenly at Tp when x > 13. 

We stopped our simulations at a; > 18, but Tp seems to continue to grow with x, and 
we see no good theoretical argument why it should not do so. Thus we conjecture that Tp 
would finally tend to infinity. 

As we said already, the average total energy was not very informative. More interesting 
was the contact energy per monomer (m) = (#(nn pair contacts)) /A^. Instead of showing it 
as a function of for various values of (x, T), we present it in fig. 8 as a function of x, at fixed 
A^ and T. We expect two clear situations, one at T > Ttripic and the other at T << Teifi). 

• In the high-T phase, T > Ttripie, we expect that (m) decreases with x, until the freezing 
line is reached. At this points (m) should have a discontinuity for N = oo. For finite 
A^ this jump should be smeared, but it should still be rather sharp, and it should not 
show any precursor. As long as the chain is in the coil phase and not close to the 
^-line, increasing stiffness should decrease the number of contacts. 

• When T is lower than Tg{0), it is also lower than the collapse temperature for stiffness 
X > 0. In this case we are from the beginning in the collapsed phase. When T << 
Tg{0), the density is rather high for all x, and we do not expect (m) to decrease initially. 
Instead, we expect first a weak increase with x, which accelerates when we pass through 
the freezing line. This time the freezing line is however much less sharp, and we expect 
much stronger finite- A^ corrections. 

Both predictions are supported by fig. 8. Most impressive is the decrease and subsequent 
jump for large T. The behavior at low temperatures is less clear, as there are still strong 
small-A^ corrections. In particular, there is an initial decrease with x for small A^. Most 
difficult to interpret are data at intermediate T. The data for T = 3.811 cross the ^-line 
at X ~ 1. For A^ = cxo we should have an increase there with infinite slope, as the specific 
heat diverges logarithmically at Tg. But the data show a steady and systematic decrease. 
Obviously this is a finite- A^ effect which dominates completely the behavior up to extremely 
large values of A^. When we approach the folding line (m) finally increases, with an A^- 
dependence intermediate between the two previous cases. 

5 Discussion 

We have been able to map out a large region of the phase diagram for semi-stiff chains 
with self avoidance and nearest-neighbor attraction. This region contains the coil-globule 
transition (second order), the freezing transition from molten globule to a folded state (first 
order), and a triple point where these transition lines meet. We found our algorithm very 
efficient, in particular for very stiff chains where we had expected the biggest problems when 
we started this investigation. 

In our study we observed three features of the phase diagram which were not expected 
on the ground of mean field theory, although the most important one, i.e. the existence of 
a triple point, was already conjectured in 



1. The collapse temperature Tg{x) is an increasing function of the stiffness; 

2. The freezing temperature Tf{x) is an increasing function of the stiffness and it does 
not show to attain any asymptotic finite value; 

3. At a critical value of the stiffness, Tp becomes higher than Tg, and the freezing tran- 
sition happens without an intermediate globular stage. 



It is remarkable that recently a phase diagram strikingly similar to fig. 3 was found in 
a (variational approximation to a) model of random copolymers [jll|. In this model, the 



freezing temperature shows a power law dependence on the variance Ae of the monomer 
pair potentials, Tp ~ Ae", in astonishing analogy with our numerical result Tp ~ x^'^^ 
(although the latter most likely does not give the correct behavior at x — > 0). It is not clear 
whether this coincidence is fortuitous or has a deeper meaning. But in any case it makes the 
present model more interesting as a toy model for protein folding. The frozen phase of the 
present model is too ordered to be taken as a good model of the native state of a protein 
(most of the bonds are parallel), but it is conceivable that the frozen disorder of amino acid 
sequences plays a similar role as the stiffness included in the present model. 

A more direct application of the present model might be to very long semiflexible polymer 
chains such as DNA or actin. Of course such polymers do not live on lattices. Thus they can 
be deformed continuously, while only discrete deformations are possible in our model. It is 
well known that going from a continuous to a discrete system can have a big effect on phase 
transitions. For spin systems, the Mermin- Wagner theorem says that this effect is mainly 
confined to 2 dimensions, but it is not a priori obvious that the same is true in the present 
case. 

Finally, we have verified that the used algorithm, PERM, is an excellent tool for studying 
polymeric systems at very low energies where all other known methods fail. Indeed, we have 
since used it to find ground states in lattice polymer models where low-lying states have been 



given in the literature. In these simulations, results of which will be published elsewhere |]22 
we have been able to find putative ground states in all cases. In several of these cases, we 
also found states lower than these putative ground states. 
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Figure 1: i?^^ 
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/N for chains with x = 10, and for Boltzmann factors q 
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Figure 2: Partition sum differences for x = 10 and L = 64, and for the eight highest energies 
shown also in fig.l (from bottom to top) 




Figure 3: Phase diagram 
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Figure 4: Fraction fg of straight joints for a; = 3 and for different temperatures, plotted against 
chain length N. At the freezing temperature, it should jump from (1 + A/q^)~^ to a constant 1. 
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Figure 5: Same as fig.4, but for x = 10. The temperatures are the same as in fig. 1 (from bottom 
to top). The freezing transition seems to take place near Tp = 4.28 it 0.02. 
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Figure 6: Anisotropy parameter p for x = 3 and for the same temperatures as in fig.4. At the 
freezing temperature, its behavior should change from power decay to a constant ^ 1. 




Figure 7: Same as previous figure, but for x = 15. At this stiffness tlie transition from 
coil to frozen state takes place directly. From bottom to top the curves represent T = 
5.613, 5.435, 5.3625, 5.339, 5.315, 5.292, and 5.246. The estimated value of Tp is 5.335 ± 0.01. 
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Figure 8: Average number of contacts (m) per monomer as a function of the stiffness for fixed 
temperatures. Three different temperatures are considered: T = 3.332 (below Tg), dashed lines; 
T = 3.811 (between Tq{Q) and Ttiipic), solid lines; and T = 5.485 (above Ttripie), dotted lines. For 
each value of T we show data for 3 different chain lengths {N = 1000, 3000, and 5000), since 
finite size effects are very important. In particular, one sees that minima of (m) are reached at 
A?^-dependent values of x, and the effective freezing temperature strongly increases with N . 



